Random restoration projects

Summary plots

# total random sites to create
tot <- nrow(restdat)

id <- stri_rand_strings(tot, 4)
dts <- range(restdat$date)

# rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
# rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
ext <- bbox(tbpoly)

lon <- runif(10 * tot, ext[1, 1], ext[1, 2])
lat <- runif(10 * tot, ext[2, 1], ext[2, 2])
tmp <- SpatialPoints(cbind(lon, lat), 
                     proj4string = crs(tbpoly)
) %>% 
  .[tbpoly, ] %>% 
  .[sample(1:nrow(.@coords), tot, replace = F), ]

restdat_rnd <- tibble(id) %>% 
  mutate(
    date = sample(seq(dts[1], dts[2]), tot, replace = T),
    top = sample(c('hab', 'wtr'), tot, replace = T)    
  )

reststat_rnd <- tibble(id) %>% 
  mutate(
    lat = tmp$lat, 
    lon = tmp$lon
  )

resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

# map by restoration type
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)

# map by date
pbase +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)

# barplot of date counts
toplo <- restall_rnd %>% 
  group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) + 
  geom_bar() + 
  coord_flip() + 
  theme_bw() + 
  theme(
    axis.title.y = element_blank()
  ) +
  scale_y_discrete(expand = c(0, 0))

Distance to restoration sites

wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')

head(wqmtch_rnd)

## # A tibble: 6 x 5
##    stat resgrp   rnk    id     dist
##   <int>  <chr> <int> <chr>    <dbl>
## 1    47    wtr     1  cTUY 3806.092
## 2    47    wtr     2  Ryvp 6696.894
## 3    47    wtr     3  fwaq 7272.233
## 4    47    wtr     4  0lGd 7962.964
## 5    47    wtr     5  peUI 9097.059
## 6    47    wtr     6  Kf3O 9290.502

Closest

## 
# plots

# combine lat/lon for the plot
toplo <- wqmtch_rnd %>% 
  left_join(wqstat, by = 'stat') %>% 
  left_join(reststat_rnd, by = 'id') %>% 
  rename(
    `Restoration\ngroup` = resgrp,
    `Distance (dd)` = dist
  )
    
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'

# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")

# base map
pbase <- ggmap(map) +
  theme_bw() +
  theme(
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  ) +
  geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
  geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)

# closest
toplo1 <- filter(toplo, rnk %in% 1)

pbase + 
  geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest twenty percent

# closest five percent
fvper <- max(toplo$rnk) %>% 
  `*`(0.2) %>% 
  ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))

pbase + 
  geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Closest all combinations

# closest all combo
toplo3 <- toplo

pbase + 
  geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)

Summarizing effects of restoration projects

Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.

salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)

## # A tibble: 6 x 3
## # Groups:   stat [2]
##    stat     cmb     cval
##   <int>   <chr>    <dbl>
## 1     6 hab_aft 24.31485
## 2     6 hab_bef 23.64315
## 3     6 wtr_aft 24.55549
## 4     6 wtr_bef 23.84958
## 5     7 hab_aft 24.81490
## 6     7 hab_bef 24.42562

head(salchg_rnd)

## # A tibble: 6 x 4
##    stat     hab     wtr     cval
##   <int>  <fctr>  <fctr>    <dbl>
## 1     6 hab_aft wtr_aft 24.43517
## 2     6 hab_aft wtr_bef 24.08221
## 3     6 hab_bef wtr_aft 24.09932
## 4     6 hab_bef wtr_bef 23.74636
## 5     7 hab_aft wtr_aft 24.86153
## 6     7 hab_aft wtr_bef 24.65012

Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.

wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)

## # A tibble: 4 x 5
##       hab     wtr              data       crv                    prd
##    <fctr>  <fctr>            <list>    <list>                 <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>

Discretization of salinity conditional probability distributions:

salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd

## # A tibble: 8 x 5
##       hab     wtr      qts       brk  clev
##    <fctr>  <fctr>    <dbl>     <dbl> <dbl>
## 1 hab_aft wtr_aft 26.39854 0.4506171     1
## 2 hab_aft wtr_aft 29.77876 0.8903498     2
## 3 hab_aft wtr_bef 26.19278 0.4321444     1
## 4 hab_aft wtr_bef 29.82180 0.8790707     2
## 5 hab_bef wtr_aft 26.35804 0.4415560     1
## 6 hab_bef wtr_aft 29.83229 0.8784237     2
## 7 hab_bef wtr_bef 26.15228 0.4249726     1
## 8 hab_bef wtr_bef 29.87533 0.8678202     2

A plot showing the breaks:

toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>% 
  unnest
ggplot(toplo, aes(x = cval, y = cumest)) + 
  geom_line() + 
  geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr) +
  theme_bw()

Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.

# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')

# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>% 
  group_by(hab, wtr) %>% 
  nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>% 
  rename(
    salev = cval.y, 
    cval = cval.x
  ) %>% 
  group_by(hab, wtr) %>% 
  nest %>% 
  left_join(salbrk_rnd, by = c('hab', 'wtr')) %>% 
  mutate(
    sallev = pmap(list(data, levs), function(data, levs){
      # browser()
      out <- data %>% 
        mutate(
          saval = salev,
          salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          salev = as.character(salev)
        )
      
      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')

chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')

chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>% 
  print(n = nrow(.))

## # A tibble: 24 x 6
##        hab     wtr salev       qts       brk  clev
##     <fctr>  <fctr> <chr>     <dbl>     <dbl> <dbl>
##  1 hab_aft wtr_aft    lo 12.385300 0.4647431     1
##  2 hab_aft wtr_aft    lo 17.290408 0.9175791     2
##  3 hab_aft wtr_aft    md  7.934825 0.4577565     1
##  4 hab_aft wtr_aft    md 10.001101 0.8997290     2
##  5 hab_aft wtr_aft    hi  4.118044 0.2745101     1
##  6 hab_aft wtr_aft    hi  5.023006 0.7254158     2
##  7 hab_aft wtr_bef    lo 14.069949 0.5312290     1
##  8 hab_aft wtr_bef    lo 19.673450 0.9323183     2
##  9 hab_aft wtr_bef    md  8.242412 0.3989296     1
## 10 hab_aft wtr_bef    md 10.204171 0.8325231     2
## 11 hab_aft wtr_bef    hi  4.205353 0.2998067     1
## 12 hab_aft wtr_bef    hi  4.942181 0.7312806     2
## 13 hab_bef wtr_aft    lo 11.636849 0.3631110     1
## 14 hab_bef wtr_aft    lo 15.729469 0.8169240     2
## 15 hab_bef wtr_aft    md  8.187762 0.4362466     1
## 16 hab_bef wtr_aft    md 10.192796 0.8625787     2
## 17 hab_bef wtr_aft    hi  4.049786 0.2586025     1
## 18 hab_bef wtr_aft    hi  4.789796 0.7109169     2
## 19 hab_bef wtr_bef    lo 12.461799 0.3684596     1
## 20 hab_bef wtr_bef    lo 17.357014 0.8163132     2
## 21 hab_bef wtr_bef    md  9.582194 0.4428301     1
## 22 hab_bef wtr_bef    md 12.651273 0.8677568     2
## 23 hab_bef wtr_bef    hi  4.137095 0.2927909     1
## 24 hab_bef wtr_bef    hi  4.708972 0.7124084     2

Final combinations long format:

chlbar_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  mutate(
    data = map(data, function(x){
      
      brk <- x$brk
      out <- data.frame(
        lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
      )
      
      return(out)
      
    })
  ) %>% 
  unnest %>% 
  gather('chllev', 'chlval', lo:hi) %>% 
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi')),
    chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
  )
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')

chlbar_rnd %>% 
  print(n = nrow(.))

## # A tibble: 36 x 5
##        hab     wtr  salev chllev     chlval
##     <fctr>  <fctr> <fctr> <fctr>      <dbl>
##  1 hab_aft wtr_aft     lo     lo 0.46474308
##  2 hab_aft wtr_aft     md     lo 0.45775646
##  3 hab_aft wtr_aft     hi     lo 0.27451009
##  4 hab_aft wtr_bef     lo     lo 0.53122896
##  5 hab_aft wtr_bef     md     lo 0.39892956
##  6 hab_aft wtr_bef     hi     lo 0.29980671
##  7 hab_bef wtr_aft     lo     lo 0.36311099
##  8 hab_bef wtr_aft     md     lo 0.43624662
##  9 hab_bef wtr_aft     hi     lo 0.25860246
## 10 hab_bef wtr_bef     lo     lo 0.36845964
## 11 hab_bef wtr_bef     md     lo 0.44283010
## 12 hab_bef wtr_bef     hi     lo 0.29279089
## 13 hab_aft wtr_aft     lo     md 0.45283601
## 14 hab_aft wtr_aft     md     md 0.44197255
## 15 hab_aft wtr_aft     hi     md 0.45090572
## 16 hab_aft wtr_bef     lo     md 0.40108937
## 17 hab_aft wtr_bef     md     md 0.43359351
## 18 hab_aft wtr_bef     hi     md 0.43147385
## 19 hab_bef wtr_aft     lo     md 0.45381302
## 20 hab_bef wtr_aft     md     md 0.42633209
## 21 hab_bef wtr_aft     hi     md 0.45231442
## 22 hab_bef wtr_bef     lo     md 0.44785357
## 23 hab_bef wtr_bef     md     md 0.42492674
## 24 hab_bef wtr_bef     hi     md 0.41961753
## 25 hab_aft wtr_aft     lo     hi 0.08242091
## 26 hab_aft wtr_aft     md     hi 0.10027099
## 27 hab_aft wtr_aft     hi     hi 0.27458419
## 28 hab_aft wtr_bef     lo     hi 0.06768166
## 29 hab_aft wtr_bef     md     hi 0.16747694
## 30 hab_aft wtr_bef     hi     hi 0.26871943
## 31 hab_bef wtr_aft     lo     hi 0.18307600
## 32 hab_bef wtr_aft     md     hi 0.13742130
## 33 hab_bef wtr_aft     hi     hi 0.28908311
## 34 hab_bef wtr_bef     lo     hi 0.18368680
## 35 hab_bef wtr_bef     md     hi 0.13224315
## 36 hab_bef wtr_bef     hi     hi 0.28759158

Discretesize chlorophyll data, all stations:

# discretize all chl data by breaks 
chlbrk_rnd <- chlbrk_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>% 
  group_by(hab, wtr, salev) %>% 
  nest %>% 
  full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>% 
  mutate(
    lev = pmap(list(data, levs), function(data, levs){
      
      out <- data %>% 
        mutate(
          lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
          lev = as.character(lev)
        )

      return(out)
      
    })
  ) %>% 
  dplyr::select(-data, -levs) %>% 
  unnest %>% 
  rename(
    chlev = lev, 
    chval = cval
    )

save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')

A bar plot of splits:

ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  facet_grid(hab ~ wtr) +
  theme_bw()

A plot showing the breaks:

toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>% 
  unnest %>%
  mutate(
    salev = factor(salev, levels = c('lo', 'md', 'hi'))
  )
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) + 
  geom_line() + 
  geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
  geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
  facet_grid(hab ~ wtr, scales = 'free_x') +
  theme_bw()